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Abstract 

This paper examines a discrete-time queuing system with applications to telecom- 
munications traffic. The arrival process is a particular Markov modulated pro- 
cess which belongs to the class of discrete batched Markovian arrival processes. 
The server process is a single server deterministic queue. A closed form exact 
solution is given for the expected queue length and delay. A simple system of 
equations is given for the probability of the queue exceeding a given length. 
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1. Introduction 

This paper provides a solution for the expected queue length and probability 
of a given queue length for a simple discrete-time queuing system. The queuing 
system in question processes one unit of work in one unit of time. Work arrives 
in integer units according to an arrival process with the following properties. 

• The system has two states, on and off. In an off state, no work will 
arrive. 

• If the system is in an off state then, with probability /o, in the next time 
unit the system is also an off state. 

• If the system is in an off state then, with probability /j, in the next time 
unit the system will move to an on state which will last for exactly i time 
units and then move to an off state. 

• If the system is on then a non-zero integer number of work units will arrive 
in this time period. The number of units of work which arrive is an iid 
random variable with g n as the probability that exactly n units of work 
will arrive in this time period. 
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This process can be modelled as a Markov- modulated process (MMP) which 
is completely characterised by the parameters fi and gi. This model has been 
studied by authors in several different areas, for example statistical physics (IBj ]. 
the study of dynamical systems Q and modelling telecommunications traffic 
[5], [16]. In the last two papers, the model is considered as the source of input 
to a network and hence it is natural to consider the queuing properties of such 
a model. In this paper expressions are derived for the expected queue length at 
equilibrium and the probability that the queue has a given length at equilibrium 
under certain natural restrictions (for example, that the utilisation of the system 
is less than one). It is shown that the expected queue length is a function of 
only four variables, the first and second moments of the parameters /; and gi. 
From Little's law [13|, the expected delay is proportional to the expected queue 
length. 

In section [2] the model is introduced formally and some basic properties are 
derived. The model is related to existing work in queuing theory. In section 
[3] the model is solved to get equations for the expected queue length and the 
probability of a given queue length. In section [5] the derived equations are 
compared with computer simulations. 
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2. A Markov queuing model 

Consider the discrete time process described in the introduction. The moti- 
vation behind this process is the idea that the lengths of off periods are mem- 
oryless (the probability that an on period begins is independent of the length 
of the off period so far) but on periods have lengths which are iid with an 
arbitrary distribution (within certain feasible constraints) . This arrival process 
is then used as the input to a queue which can process one unit of work (either 
queued or newly arriving) per time period. 

The process Y t , the arrival process to the system, can be modelled by an 
MMP. The process has the underlying Markov chain shown in figure [TJ At any 
time the chain is in an on state (all of the non-zero states) then a non-zero 
number of arrivals will occur with a given probability. 

The n + 1 state Markov chain from figure [T] has the transition matrix 
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Figure 1: The Markov chain for an MMP traffic model. 

Let {X t : t S N} be the states of a discrete-time, homogeneous Markov 
chain at time t. This chain is shown in figure [T] and the parameters fi define 
the transition probabilities of the chain. It can be seen from the structure of 
the chain that /, is the probability that, given the chain is currently in state 
zero, the next state will be i. This could also be thought as the probability 
that, given the chain is currently in an off state, the next state will begin an 
on period of length exactly i (/q is the probability the next state will be an off 
state again). 

Let gi : i £ N be the probability that exactly i units of work arrive when the 
chain is in an on (non-zero) state (X t > 0). Let {Y t : t s N} be a series derived 
from Xt by the rule, that if Xt = then Y t — and if Xt > then for all 
i G {1, 2, . . . , m} (where m is the maximum possible value of Yt) P [Y t = i] = gi 
(for an on state, the possibility that Y t = is excluded, that is, go = 0). 

The server model used is a deterministic process where exactly one unit of 
work is processed in one time unit. 

2.1. Related work 

The arrival process belongs to the processes known as discrete batch Marko- 
vian arrival process (D-BMAP) and the queuing system is a subset of D-BMAP /D/l 
queues. The BMAP itself was introduced by Lucantoni Both the D-BMAP 
and the BMAP have previously been used as models of telecommunications traf- 
fic i,0Ji- Details of the BMAP can be found in most modern books on queuing 
theory [J, Chapter 12] and only a brief outline is given here. It is also interesting 
to note that the arrival model considered here is very similar to the batch re- 
newal process studied by Fretwell and Kouvatsos in the context of internet 
traffic. In fact the system they study is the one in figure Q] with on and off 
reversed. 

The structure of a generic D-BMAP is as follows. Let P be the transition 
matrix for a discrete time Markov chain with state space E = Nx 0, . . . ,n. 
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Each state of the chain is a pair (J, i) where j is the level (the number of 
arrivals generated in that state) and i is the phase. The transition matrix has 
the structure 

D Di D 2 D 3 
D Di D 2 
D Di 



where D^ is, itself an (n + 1) x (n + 1) matrix. The (j, k)th entry in D^ is the 
probability of a phase transition from phase j to k given level i. 

The arrival process described in the previous section could be described as 
a D-BMAP where the phase is the state of the chain in Figure Q] and the levels 
simply represent the various on states. The D are the (n+1) X (n+ 1) matrices 
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for i 6 {1, 2, ... , to}. 

Blondia and Cassals 0] provide a method for solving the D-BMAP/D/1/K 
queuing model which gives a solution in terms of a recursive series of matrix 
equations. The complexity of calculation is given by the author as K 2 M 3 where 
K is the buffer size (the possible number of items in the queue) and M is the 
number of phases (n + 1 in the system described here). 

In contrast the system studied here only works for infinite buffers and gives 
an answer for the expected queue length in closed form. It gives an answer 
for the probability of the queue having a given length as a recursive system of 
equations each with n terms. 



2.2. Basic properties of the system 

It is useful to define / = J27 =1 ifi and f 2 — X)"=i* 2 /«- Similarly, define 
9 = Sl^i *5» an d a ^ so 9 2 = Eti i 2 9i- Since g is the mean number of arrivals in 
an on state then this should be finite (otherwise the system will have an infinite 
mean arrival rate). This will be the case in all systems with to finite. 

Let 7Tj be the equilibrium probability of state i. This exists when the chain 
is ergodic. It can be easily shown that the finite chain is ergodic if fo G (0, 1) 
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and, if this is the case, the equilibrium probability of state i, 7r^ is given by 



and 

" 1 

and rearranging gives 

It is also of interest to consider the version of this chain where n — > oo. That 
is for any JVgN there exists i > N such that > 0. The mean return time 
to the state zero is given by 1 + / which is finite if / is finite. Therefore the 
infinite chain is ergodic if /o > and / is finite. The infinite chain is useful in 
studying long-range dependence P, Obviously for finite chains / is finite. 

The output Y t is then used as input to a queuing system. Let Q t represent 
the number queuing at time t and Y t represent the number of arrivals to the 
system during [t, t+1). One queued item is removed from the queue every time 
unit if there is any work in the queue or if any arrives (if Qt + Yt > 0) . Therefore 
the queuing system is as follows 

Qt+i = [Qt +Y t - 1] + , (1) 

where [X] + means max(0,X). 

The quantity 1 — ir — //(l + /) represents the proportion of the time the 
system is in an on state and therefore g(l — t:q) is the mean arrival rate A. Since 
the system can output one unit of work per time period the queue utilisation 
(proportion of time the queue is non-empty) p = A. Both are given by, 



gf_ 



(2) 



In order that the queue does not grow forever it is obviously a necessary condi- 
tion that p < 1. 

To summarise, the model is specified by the parameters /j and gi. The 
requirements on these parameters are that /o > (which guarantees the under- 
lying MC is aperiodic), that / is finite and that gf/(l + /) < 1 (which in turn 
is a requirement that g is finite) . The first two requirements together guarantee 
that the underlying MC is ergodic the third requirement ensures that the queue 
utilisation is less than one. 

2.3. Notation 

The following notation is used in this paper and is gathered here for conve- 
nience. 
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• Y t — the number of arrivals to the system during [t,t + 1). 

• X t — the state of the underlying Markov chain during [t,t + 1). 

• Qt — the number queuing at time t. 

• X,Y,Q — the above quantities as random variables at some time when 
the system is in equilibrium. 

• TTi — the equilibrium probability of the ith state of the chain. 

• P — the transition matrix for X t . 

• fi — the transition probabilities in P. 

• n — the highest possible value of X t (the highest numbered state in P). 

• /> f 2 — the first and second moments of fi, Y17=i ifi an d T^i=i i 2 fi- 

• Qi{z) — the generating function J^^Lo ^ [Q ^ Qi X ~ i] z q . 

• Q(z) — the n + 1 column vector [Qo(z), Qi(z), . . . , Q n (z)] T . 

• 9i — the probability that an amount of work i arrives in the next time 
unit if the system is on, P [Yt — i\Xt > 0]. 

• 77i — the maximum possible value of Y t (the largest number of units of 
work which may arrive in unit time) . 

• 5i <? 2 — the first and second moments of <?;, anc ^ ^2T=\ i 2 9i- 

• g(z) — the generating function Y^iLi 9i zl ■ 

• bo — the boundary condition P [Q = 0\X = 0]. 

• p — the queue utilisation. 

• A — the mean arrival rate. 

3. Solving the queuing model 

Let Q(z) = {Qo(z), Qi(z), . . . , Q n {z)] T be a column vector of the generating 
function for the queue in each state of the chain. That is, 

oo 

Q l (z)=Y,V[Q = q,X =i] z q . 

q=0 

Consider a general MMP arrival process. Let Ai(z) be the generating func- 
tion for the number of arrivals if the underlying chain is in state i. Let B = 
[bo, bi, . . . , b n ] T be the (n + 1) column vector of boundary conditions, bi — 
P [Y = 0, Q = 0\X = i]. Following Li it can be shown that the queuing 



G 



system of equation ([TJ) using a general MMP with transition matrix P as input 
implies 

Q(z) = (z-l)[zI-P T G(z))]- 1 P T B, (3) 

where G(z) = diag(A (», A\{z), . . .,A n (z)). 

For the MMP in question B and G(z) have much simpler forms. In the 
off state the generating function for arrivals is simply 1 (no arrivals occur with 
probability one) and in the on state the generating function is g(z). There- 
fore, G(z) = diag(l, g{z), g(z), . . . , g(z)). Since in the on state, the system 
always has at least one arrival and in the off state the system always has 
no arrivals then B is the (n + 1) column vector, B = [bo, 0, 0, . . . , 0] T where 
b Q = P [Y = 0, Q = 0\X = 0] = V[Q = 0\X = 0]. It is these simplifications 
which make this system soluble. 

For the specific MMP being studied ((3]) can be rewritten, 



Q(z) = («-1)(A- 
where A is the (n + 1) x (n + 1) matrix 
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and A' is the rank one (n+1) x (n+ 1) matrix which can be written as A' = uv 
where u = [f , /i,/ 2 , . • . ,/„] T and v = [-1,0,0, . . .]. 

The matrix A can be inverted giving the (n + 1) x (n+1) matrix 



A" 1 



(g(z)/zf (g(z)/zY (g(z)/zf {g{z)/zf 

(g{z)/zf (g(z)/zY (g{z)/zf 

{g(z)/zf (g(z)/zY 

(g(z)/z)° 



Note that since .go = then g{z)/z) = Y^iLi 9i z% I z = Y^Li 9i z% 1 ■ If z G [0, 1] 
then £™i S^ 1 " 1 < YTiUdi = L Therefore g{z)/z G [0,1] if z G [0,1] Hence 



[g(z)/z) n remains bounded as n 
for example, Q) states that, 



oo. The Sherman-Morrison formula (see, 



(A + AT^A-'- fY," 1 
1 + vA-^u 



Now, 



1 + vA-\ = 1-1/ z^2(g(z)/zyfi. 

i=Q 
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Define 



Then 



whence, 



(A + A 



j=i 

[1 - a Q (z)]I - A^uv 



1 - a (z) 



A"\ 



Q(z) = (z-l) 
Multiplying the matrices gives 



[l-o (z)]I-A ^ A -i p T B 
l-a (z) 



Qt{z) 



z - 1 



-b ai(z). 



1 - a Q (z) 

From the definition of Qi(z), E \z®\ — J27=o Qi( z ) an d therefore 

E[z]=6 ° -EL/^ = ^P () 

where N(z) = (z-l) Elo fi YtjM*)/*)' and D ^ z ) = z ~ E?=o fM*)/*?- 

3.1. Calculating bo 

Now it is necessary to calculate &o- Note that lim^i z® = 1, whence 

z->i iV(z) 

Since, lim z _i -D(z) = lim z _»i N(z) = 0, by L'HopitaPs rule, 

bo = Inn ^7^-- 

Z->1 iV'(z) 



But 



and 



C'(z) = l-(.9(z)/z)'^z/,( 5 (z)/zr 1 , 



i=l 



{g(z)/z)' = z 2 - ft*') ■ 

i=l 

Hence lim z ^i(g(z) / z)' = ~g — 1, which implies, 

n 

lim =l-(3-l)Vi/ i = l + (l- #7. 

2— »1 z — ' 



(5) 



i=l 
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Similarly, 

n i n i 

N '( z ) = E /.EW 2 )/^ + (* - !) 'EteWM'fi E M*)/*)*' 1 - 

i=0 j=0 i=0 j=0 

Providing the sum at the right hand side remains finite (which it will for all 
finite n) then the factor of [z — 1) will cancel this term as z — > 1. This gives 



Urn N'{z) = V(i + = 1 + /. 

z— >1 * ' 

i=0 

Finally, therefore, 

b Q = l-^L = l-p. (6) 

<?.#. Queuing results 

The next stage is to get a function for the expectation of the queue size. 

dE \z Q ] °° °° 

Hm L — L = n m q p [Q = q ] z 9-i = V qP [Q = q] = E [Q] . 

z— >i az z— ' ^ — ' 

q=l q=l 

Since &o is constant, from (HI), 

d(N(z)/D(z)) N'(z)D(z)~N(z)D'(z) 
E [Q] = &0 Inn = b hm ^ . 

Similarly, lim z ^i D(z) — and, from L'Hopital's rule, 

EtQl^olim^^-^^ 
Ml 2D(z)D'(z) 

and, since lim z ^i N(z)/D(z) = l/&o, 

It is now necessary to find expressions for N"(z) and D"(z). Firstly, 
AT"(z) = POK*)/*)' + (* - l)Cff(*)/*)"] H=i /i Ej=i M*)/*)*- 1 

+(z i)(. 9 (z)/z)' 2 EILi A£* - l)(.g(z)/zp- 2 . 

Hence, if all the sums remain finite (which they will if n is finite), 

n 

hm N"(z) = [2(g - 1)] £ i(i + l)/ 4 /2 = (9-l)(7 5 + 7)- (8) 
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Similarly 



D"{z)=-(g{z)/z)"Y d if i {g{z)/z) i - 1 
i=l 

n 

-(gizyzy^iii-lMgiz)/; 



ii-2 



and 



lim( 5 (z)/z)" = lim z~ 3 V(i 2 - 3z + 2)g i z i = g 2 -3g + 2, 



therefore, 

Substituting 



' 5 • f <r 
and © into gives 



272 



lim !>"(*) = [5 + r - .9 2 - 1]/ - (9 ~ Iff 

z—>l 



E[Q] 



5(5-1) 



./• 2 



7 2 



+ /(! + /) 



2(l + f)[l + f-fg] 



(9) 



(10) 



—2 



Note that 5 > 1 , / 2 > /' and g 2 > 7? 2 by their respective definitions and 1 + 
f > fg for a system with utilisation less than one by equation @ . All bracketed 
terms in the numerator are therefore positive or zero and the denominator is 
strictly positive. Note, however that f 2 and g 2 are only guaranteed finite for 
systems with finite n and m respectively. That is to say that some systems 
with a mean arrival rate less than one and an utilisation less than one will still 
have an no finite value for the expected queue length. An example of such a 
system are the systems with the fi parameters given in 3, B[ which both have 
f 2 as a non-convergent series. Such systems are of interested to those studying 
long-range dependence and heavy-tailed distributions. 



Interpreting f 2 — f and g 2 
T0|) could also be written as 

E [Q] = 



g as the variance of / and g respectively then 



-l)var(/) + /(l + /)var( ff ) 



2(l + /)2(l-p) 
From |J2J|, the expected delay is given from Little's law 



E [T] 



E [Q] _ g{g - l)var (/) + /(l + /)var (g) 



A 



2(1 + /)V(1-P) 



An implication of these equations is that, assuming the mean traffic level 
is fixed (that is / and ~g cannot be changed) then the queueing delay of the 
system would be minimised if the variance in the lengths of the on periods was 
minimised and the variance of the amount of traffic arriving in an on period 
was minimised. This has an interesting analogy to the well-known Pollaczek- 
Khinchin result Q that for an M/G/l queue the waiting time is proportional 
to the variance in the service time. 
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3. 3. Finding the queue distribution function 

In order to be able to ask questions about, for example, buffer overflow 
probabilities, it would be useful to be able to ask questions about the probability 
of a given queue size (P [Q = i]) or the probability that the queue is more than 
a given size (f[Q > i]). 

The probability that the queue is zero is given by 

P [Q = 0] = Urn E [z Q ] , 

and, more generally, the probability that the queue length is q can be found by 
differentiating q times and taking the limit as z —* 0. 

d q E \z Q ] 
PQ = i=lim 1 J . 

z^o q\diz 

This can be solved computationally by repeated symbolic differentiation. How- 
ever, this is computationally intensive and the algorithm is numerically unstable. 
Another approach is to produce a recursive formula for the coefficient of z l in 
E [z®] . This can be done using standard techniques for formal power series 
from, for example, Knuth 

By definition, the coefficient of z 1 in E [z®~\ is P [Q = i] . Let iVj and Di be 
the coefficients of z % in N(z) and D(z). Since (j4|) is true for all z, therefore 
standard techniques for division of power series give 

fe 

£>[Q = i] D k -i = b N k , 

i=0 



which rearranges to 



P \Q = k] = — 



k-l 



N k b - p [Q = *] D 



k—i 



(ii) 



This recursive formula expresses P [Q = k] in terms of bo which can be evaluated 
with ©, coefficients N t and Di and P [Q = j] for j < k. Therefore, the P [Q = k] 
can be calculated in turn beginning with P [Q = 0] which is given by 

P[Q = 0] = M£. 

Uo 

The coefficients Ni and Di can be easily calculated. Let Gij be the coefficient 
of z % in (g(z)/ z)i . The coefficients above can be expressed as 

n 

D i = 5i- 1 -^f j G i j (12) 

3=0 

where Si is the Kronecker delta function (#j_i = 1 if i = 1 and otherwise) and 

N , I - E,=o T,k =J fk * = o , . 

1 lv; .,:<;, ,, a,., y.1 ,.a ■ ■•». [ 1 
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The coefficients Gi.j can be caiculated by another recurrence relation. Since 



9{z)/z = 9i+iz\ 
and as G^o = &i> this gives the recurrence relation 

i 

Gi,j+1 = / ] Gk,j9i+l-k- 
fc=0 



3.4. A simpler model — the constant batch size model 

A simplification occurs when the batch size is fixed. Assume that, work must 
arrive in units of exactly r where r > 1 (r = 1 is the uninteresting system where 
no queue ever forms). This means that g(z)/z = z r ~ 1 and (g(z)/zy = z^ r ~ 1 '. 



In turn this gives G 



1,3 



where here, and throughout this section, 5 is 



the Kronecker delta function. Obviously g = r and g 2 

m _ 

r(r-l)(/2 
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hence from 



E [Q] = 



2(l + /)(l + /-r/) 

bo/ fo and for k > 0, 



k] = T 
Jo 



fc-i 



It can also be shown that P [Q = 0] 

— 6o^(fc-l)/(r-l)-L(fc-l)/(r-l)J 
+ &0<5fc/(V-l)-|fc/(r-l)J 



Q = k - 1] - 5 i/(r-i)-u/(r-i)j/j p [<5 = i] 



j=(fc-l)/(r-l) 



E 

j=fc/(r-l) 



where [xj is the floor function and with the notational conveniences that fj = 
for j > n and that ^™ /j = for j > n. Note that the delta functions here are 
simply testing if a given expression, for example k/(r — 1), is integer. From ([B]), 



b 



1 + / ~ rf 
1 + f 



4. Simulation tests 

The system needs to be tested against simulation to see if it can be practically 
used. While the equations from the previous sections are correct, they are not 
useful if the numerical stability of the recursive system of equations is poor. 
For the expected queue length calculations this is not an issue but it is for the 
probability of higher queue lengths since it is likely that P [Q = i] will become 
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extremely small as i becomes large. This, in turn, will make the calculation of 
(fTTj) problematic as i becomes large. Calculating the potential effects of errors 
in the system of equations given by (fTTj) , (TT2"|) and (Tj"3"| is non-trivial. The 
answers presented here are tested against simulation methods and appear valid 
for smaller queue sizes but become obviously incorrect (negative probabilities for 
example) for larger queue sizes. The simulations here were done in python. The 
same calculations have been tried using arbitrary precision arithmetic libraries. 
This enabled slightly larger queue sizes to be calculated and give reasonable 
answers but at the expense of greatly increased run time. 

The method is simply to replicate the Markov chain and queuing system 
described earlier and simulate it. This is an exact simulation of the queuing 
system described in the paper. The simulation can then run for a set number of 
iterations and the queue measured at each point to sample E [Q] or get a sample 
of the probabilities P [Q = i] by measuring the proportion of the iterations where 
the queue has the value i in simulation. Tables [JJ or [2] show the parameter sets 
for two different simulation scenarios with the first representing a lightly loaded 
system and the second representing much heavier loading. In both scenarios 
E [Q] matches well between theory and experiment as expected and no results 
are presented here. 



i 


12 3 


fi 
9i 


0.8 0.1 0.05 0.05 
0.4 0.4 0.2 



Table 1: Parameter set 1 for simulation. 



i 





1 


2 


3 


4 


fi 


0.6 


0.2 


0.1 


0.05 


0.05 


9i 





0.2 


0.6 


0.1 


0.1 



Table 2: Parameter set 2 for simulation. 

Figures [2] and [3] show the theory from section [3~3l plotted against ten simula- 
tion runs for each parameter set, each run being 10 s iterations. The simulation 
results are shown as a mean for the ten simulation runs and upper and lower 95% 
confidence intervals. The plots are on a logscale and therefore, obviously, values 
which are zero or negative do not show up. This happens with lower confidence 
intervals and when there are rounding errors in the theoretical calculation. 

For parameter set one, rounding errors make the theoretical calculation obvi- 
ously unreliable (this is obvious when the numbers become negative) at around 
i = 40 (where P [Q = i] is around 1 x 10 -16 ) although it is likely that some figures 
before this are rendered inaccurate due to rounding — the negative numbers are 
omitted from the plot because of the logscale, all others are included. In fact 
Note that because only 10 8 iterations were simulated the lowest sample prob- 
ability that can be assigned in simulation is 1 x 10~ 8 . For parameter set two, 
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1 




1e-18 1 1 1 1 1 — 1 1 1 1 — 1 1 1 1 1 

1 10 100 



Figure 2: Theory versus simulation for parameter set 1. 

it can be seen from figure [3] that the calculation remains reliable until above 
i = 110 when, again, it becomes obviously inaccurate due to rounding at a 
probability again around 1 x 10~ 16 . As has been mentioned, by using arbitrary 
precision arithmetic libraries, the numerical stability can be increased slightly 
but at the expense of greatly increased run times. 

5. Conclusions 

The arrival system described is quite general and could be useful in any sys- 
tem when work arrives at discrete times in discrete batches. The solutions given 
provide mean queue lengths and delays for the given arrival process. In addition 
a system of equations has been given which can calculate the probability that 
the queue has a given length from the system parameters and the probabilities 
of smaller lengths. The numerical stability of the recursive system of equations 
giving the probability distribution has been assessed via simulation. 
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